classdef lpcvd_2DSurface
    properties
        x = []; % Unit nm
        y = []; % Unit nm
        h = []; % Unit nm
        site_width = 0.2;
        site_width_unit = 'nm';
        layer_thickness = 0.2;
        layer_thickness_unit = 'nm';
    end
   
    properties(Dependent = true)
        r2 = 0;
        m2 = 0;
        autocor = 0;
    end
    
    methods
        function surface = lpcvd_2DSurface(filename)
            if(nargin==1)
                surface.h = dlmread(filename);
                surface.x = surface.site_width*(0:(size(surface.h,1)-1))';
                surface.y = surface.site_width*(0:(size(surface.h,2)-1))';
            end
        end
        
        function h = mesh(surface)
            h = mesh(surface.x,surface.y,(surface.h)');
            set(gca,'FontSize',18);
            xlabel('x (nm)');
            ylabel('y (nm)');
            zlabel('z (nm)');
            title(sprintf('r^2 = %.3e (nm^2), m^2 = %.3e',surface.r2,surface.m2));
        end
        
        function obj = set.site_width(obj,width)
            site_width_old = obj.site_width;
            obj.x = obj.x/site_width_old*width;
            obj.y = obj.y/site_width_old*width;
            obj.site_width = width;
        end
        
        function r2 = get.r2(obj)
            h_avg = mean(mean(obj.h));
            h_dev = obj.h-h_avg;
            r2 = sum(sum(h_dev.^2))/(length(obj.x)*length(obj.y));
            
            %             for i = 1:length(obj.x)
            %                 for j = 1:length(obj.y)
            %                     r2 = r2+(h(i,j)-h_avg)^2;
            %                 end
            %             end
            %             r = sqrt(r2/(length(x)*length(y)));
        end
        
        function m2 = get.m2(obj)
            s = [diff(obj.h,1,2),obj.h(:,1)-obj.h(:,end)]./(obj.x(2)-obj.x(1));
            nx = length(obj.x);
            ny = length(obj.y);
            m2 = sum(sum(s.^2))/(nx*ny);
%             for i =1:nx
%                 for j = 1:ny
%                     m2 = m2+s(i,j)^2;
%                 end
%             end
%             m = m2/(nx*ny);
        end
        
        function s2 = aggregate(s1,A)
            s2 = lpcvd_2DSurface();
            
            s2.x = s1.x(1:A:end)/A;
            s2.y = s1.y(1:A:end)/A;
            s2.site_width = s1.site_width*A;
            s2.site_width_unit = s1.site_width_unit;
            s2.layer_thickness = s1.layer_thickness;
            s2.layer_thickness_unit = s1.layer_thickness_unit;
            
            dim_s2 = length(s2.x);
            s2.h = zeros(dim_s2,dim_s2);
            for i = 1:dim_s2
                for j = 1:dim_s2
                    s2.h(i,j) = mean(mean(s1.h((i-1)*A+1:i*A,(j-1)*A+1:j*A)));
                end
            end
        end
        
        function s2 = resample(s1,A)
            s2 = lpcvd_2DSurface();
            
            s2.x = s1.x(1:A:end)/A;
            s2.y = s1.y(1:A:end)/A;
            s2.site_width = s1.site_width*A;
            s2.site_width_unit = s1.site_width_unit;
            s2.layer_thickness = s1.layer_thickness;
            s2.layer_thickness_unit = s1.layer_thickness_unit;
            
            dim_s2 = length(s2.x);
            s2.h = zeros(dim_s2,dim_s2);
            for i = 1:dim_s2
                for j = 1:dim_s2
                    s2.h(i,j) = s1.h((i-1)*A+1,(j-1)*A+1);
                end
            end            
        end

        function sp = calcSpectrum(obj,mode)
            sp = spectrum2D(mode);
               
            x = obj.x;
            y = obj.y;
            L = max(x);
            l = length(x);
            assert(2*mode<l,'Site = %d, no more than %d modes are available',mode,floor(l/2));
            
            temp = L^2/l^2;
            for m = 0:mode
                for n = 0:mode
                    phi1mn = fun_phi1mn(m,n,x,y,L);
                    phi2mn = fun_phi2mn(m,n,x,y,L);
                    phi3mn = fun_phi3mn(m,n,x,y,L);
                    phi4mn = fun_phi4mn(m,n,x,y,L);
                    sp.z1(m+1,n+1) = sum(sum(obj.h.*phi1mn))*temp;
                    sp.z2(m+1,n+1) = sum(sum(obj.h.*phi2mn))*temp;
                    sp.z3(m+1,n+1) = sum(sum(obj.h.*phi3mn))*temp;
                    sp.z4(m+1,n+1) = sum(sum(obj.h.*phi4mn))*temp;
                end
            end
            
            function phi1 = fun_phi1mn(m,n,x,y,L)
                phi1 = (2/L)*sin(2*m*pi*x/L)*(sin(2*n*pi*y/L))';
            end
            
            function phi2 = fun_phi2mn(m,n,x,y,L)
                if (m== 0 || m == l/2) && (n==0||n==l/2)
                    coeff = 1/L;
                elseif m~=0 && m~=l/2 && n~=0 && n~=l/2
                    coeff = 2/L;
                else
                    coeff = sqrt(2)/L;
                end
                phi2 = coeff*cos(2*m*pi*x/L)*(cos(2*n*pi*y/L))';
            end
            
            function phi3 = fun_phi3mn(m,n,x,y,L)
                if n == 0 || n == l/2
                    coeff = sqrt(2)/L;
                else
                    coeff = 2/L;
                end
                phi3 = coeff*sin(2*m*pi*x/L)*(cos(2*n*pi*y/L))';
            end
            
            function phi4 = fun_phi4mn(m,n,x,y,L)
                if m ==0 || m == L/2
                    coeff = sqrt(2)/L;
                else
                    coeff = 2/L;
                end
                phi4 = coeff*cos(2*m*pi*x/L)*(sin(2*n*pi*y/L))';
            end
            
        end
       
    end
end

